# -*- coding: utf-8 -*-
"""
/**************************************************************************
        SDATool - A QGIS plugin to do spatial data analysis 
                                using R Statistical package
                             -------------------
    begin                : 2008-01-15
    copyright            : (C) 2008 by Volkan Kepoglu
    email                : vkepoglu at gmail.com
 *************************************************************************/
/**************************************************************************
 *   This program is free software; you can redistribute it and/or modify *
 *   it under the terms of the GNU General Public License as published by *
 *   the Free Software Foundation; either version 2 of the License, or    *
 *   (at your option) any later version.                                  *
 *************************************************************************/
"""
# R options
import sda4ppRoptions
# write to logFile.log
import wl

def doSpatialAnalysis (r,index):
  dt="diggleRpart-doSpatialAnalysis"
  wl.wl("index: " +str(index),dt)
  # to data frame
  r("v_sp <- as(v_spdf,\"SpatialPoints\")")
  r("Fire98 <- v_sp@coords")
  r("v2_sp <- as(v2_spdf,\"SpatialPoints\")")
  r("Nofire98 <- v2_sp@coords")
  # new library
  r("library(splancs)")
  # common funcs
  r("xcoord<-c(Fire98[,1],Nofire98[,1])")
  r("ycoord<-c(Fire98[,2],Nofire98[,2])")
  r("Nofire0 <- as.matrix(rbind(Nofire98))")
  r("Nofire<-cbind(Nofire0[,1]-min(xcoord),Nofire0[,2]-min(ycoord))")

  r("fire98<-cbind(Fire98[,1]-min(xcoord),Fire98[,2]-min(ycoord))")
  r("fire.a <- as.matrix(rbind(fire98,Nofire),ncol=2)")
  r("labs<-rep(c(1,0),c(nrow(fire98),nrow(Nofire)))")
  r("fire.bb <- bboxx(bbox(as.points(fire.a)))")
  r("t <- seq(from = 1, to=100, by=1)")

  r("poly<-fire.bb")

  r("m<-mse2d(fire.a,poly,100,1000)")
  r("bw<-m$h[which.min(m$mse)]")

  r("pts1<-as.points(fire98[,1],fire98[,2])")
  r("pts0<-as.points(Nofire[,1],Nofire[,2])")

  # plot
  if index["plotOpt"] == "diggle":
    r("k11h <- khat(fire98, fire.bb, t)")
    r("k22h <- khat(Nofire, fire.bb, t)")

    r("D12hatRL<-vector()")
    r("set.seed(1)")
    r("for (i in 1:100){relab<-sample(labs);Fire<-fire.a[relab==1,];NoFire<-fire.a[relab==0,];k11 <- khat(Fire, fire.bb, t);k22 <- khat(NoFire, fire.bb, t);D12hatRL<-rbind(D12hatRL,k11-k22)}")

    r("SE.est<-apply(D12hatRL,2,FUN=sd)")

    r("Ylim<-range(c(2*SE.est,-2*SE.est,k11h-k22h))")

    rCode = ["par(mfrow=c(1,1),cex.axis=1.5,cex.lab=2,cex.main=2,mar=c(5,6,4,2))",
             "plot(t, k11h-k22h, type=\"l\", xlim=c(0,100), ylim=Ylim, main=\"Diggle D Function\", xlab=\"s(m)\", ylab=expression(hat(D)(s)),lwd=2)",
             "abline(h=0,lty=3,lwd=3)",
             "lines(t, 2*SE.est, type=\"l\",lty=2,col=\"red\",lwd=4)",
             "lines(t, -2*SE.est, type=\"l\", lty=2,col=\"red\",lwd=4)"]
    sda4ppRoptions.plotGraphics(r,rCode)
  elif index["plotOpt"] == "bandwidth":
    rCode = ["plot(m$h,m$mse,type=\"l\",main=paste(\"Bandwidth = \", bw))"]
    sda4ppRoptions.plotGraphics(r,rCode)
  elif index["plotOpt"] == "case":
    r("ker1<-kernel2d(pts1,poly,bw,250,250)")
    rCode = ["image(ker1, main=\"Intensity of Case Points\",col=heat.colors(12))"]
    sda4ppRoptions.plotGraphics(r,rCode)
  elif index["plotOpt"] == "background":
    r("ker0<-kernel2d(pts0,poly,bw,250,250)")
    rCode = ["image(ker0, main=\"Intensity of Background Points\",col =heat.colors(12))"]
    sda4ppRoptions.plotGraphics(r,rCode)
  elif index["plotOpt"] == "ratio":
    r("kr<-kernrat(pts1,pts0,poly,h1=100,h2=100,250,250)")
    rCode = ["image(kr, asp=1, main=\"Ratio of Case to Background Intensity\" )"]
    sda4ppRoptions.plotGraphics(r,rCode)

